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Abstract. In this paper, we derive the minimum-energy periodic control that 
entrains an ensemble of structurally similar neural oscillators to a desired frequency. 
The state space representation of a nominal oscillator is reduced to a phase model by 
computing its limit cycle and phase response curve, from which the optimal control 
is derived by using formal averaging and the calculus of variations. We focus on the 
case of a 1:1 entrainment ratio, and introduce a numerical method for approximating 
the optimal controls. The method is applied to asymptotically control the spiking 
frequency of neural oscillators modeled using the Hodgkin-Huxley equations. This 
illustrates the optimality of entrainment controls derived using phase models when 
applied to the original state space system, which is a crucial requirement for using 
phase models in control synthesis for practical applications. The results of this work 
can be used to design low energy signals for deep brain stimulation therapies for 
neuropathologies, and can be generalized for optimal frequency control of large-scale 
complex oscillating systems with parameter uncertainty. 
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1. Introduction 

The synchronization of oscillating systems is an important and extensively studied 
scientific concept with numerous engineering applications |T]. Examples include 
the oscillation of neurons [2], sleep cycles and other pacemakers in biology [31 
m |5], semiconductor lasers in physics [6], and vibrating systems in mechanical 
engineering [7]. Among many well-studied synchronization phenomena, the asymptotic 
synchronization of an oscillator to a periodic control signal, called entrainment, is of 
fundamental scientific and engineering importance [H |9] . The intrinsic occurrence and 
extrinsic imposition of entrainment in networked oscillators is of particular interest in 
neuroscience [lOl [11]. It has been investigated as an important mechanism underlying 
the coordinated resetting of neural subpopulations, which is required for demand- 
controlled deep brain stimulation (DBS) for clinical treatment of Parkinson's disease [T2] . 
In clinical DBS, entrainment is accomplished using a sequence of pulses with tunable 
amplitude, duration, and frequency, after which another effectively desynchronizing 
stimulation is applied. Because low energy consumption increases implant battery life, 
decreases tissue damage, and minimizes side effects and risks [13], the development of 
minimum energy stimuli for the entrainment of neural oscillator populations is critical 
to improving the clinical outcomes of DBS. 

The entrainment, and hence frequency control, of an oscillating system can be 
examined by considering its phase response curve (PRC) [HI [15], which quantifies 
the shift in asymptotic phase due to an infinitesimal perturbation in the state. The 
classic phase coordinate transformation [16] for studying nonlinear oscillators was used 
together with formal averaging [17] to develop a model of synchronization in coupled 
chemical oscillations [H]. Since then, phase models have become indispensable in 
physics, chemistry, and biology for studying oscillating systems where the full state- 
space model is complicated or even unknown, but where the phase can be estimated 
from partial state observations, and the PRC can be approximated experimentally [19] . 
They have been successfully applied to investigate many synchronization phenomena 
[20] . focusing on synchronization emerging in networks of interacting oscillators and on 
the response of large collections of oscillators to periodic external stimuli [211 122] • Such 
models have long been of interest to neuroscientists [231 [21] , who have been motivated 
by the prospect of using dynamical systems theory to improve the effectiveness of DBS 
as a clinical therapy for epilepsy and Parkinson's disease [251 ES]. Several studies have 
concentrated on the use of phase models in order to attain desired design objectives 
for electrochemical [271 [28] and neural [291 [30] systems, including recent work that 
approaches the use of phase models in neuroscience from a control theoretic perspective 
[3T| [32] . The control of neural spiking using minimum energy controls with constrained 
amplitude and charge balancing has also recently examined [331 El] • These studies have 
demonstrated that phase-model reduction provides a practical approach to synthesizing 
near-optimal controls that achieve design goals for oscillating neural systems. 

Much of the work on the control of neural oscillators is based on the assumption that 
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each neuron behaves according to pre-defined underlying dynamics, such as the Hodgkin- 
Huxley equations [35], which constitute a widely studied model of action potential 
propagation in a squid giant axon. However, in practical applications of neural control 
and engineering, the systems in question are collections of biological neurons that exhibit 
variation in parameters that characterize the system dynamics, specifically the frequency 
of oscillation and sensitivity to external stimuli. Although such a system consists of a 
finite collection of subsystems, it contains so many unobservable elements, each with 
parameter uncertainty, that its collective dynamics are most practically modeled by 
indexing the subsystems by a parameter varying on a continuum. The control of such 
neural systems therefore lies within an emerging and challenging area in mathematical 
control theory called ensemble control, which encompasses a class of problems involving 
the guidance of an uncountably infinite collection of structurally identical dynamical 
systems with parameter variation by applying a common open-loop control |36]. In 
the context of phase model reduction, the appropriate indexing parameter for such a 
collection of oscillating systems is natural frequency. Therefore, a practical approach 
to the optimal design of inputs that entrain a collection of neurons is to first consider 
a family of phase models with common nominal PRC and frequency varying over a 
specified interval. Optimal waveforms that entrain a collection of phase oscillators with 
the greatest range of frequencies by weak periodic forcing have been characterized for 
certain oscillating chemical systems |37], and this approach has been extended to a 
method for optimal entrainment of oscillating systems with arbitrary PRC [38] . 

In this paper, we develop a method for engineering weak, periodic signals that 
entrain ensembles of structurally similar uncoupled oscillators with variation in system 
parameters to a desired target frequency without the use of state feedback. In addition, 
we present an efficient numerical method for approximating optimal waveforms by 
minimizing over a compact interval a polynomial whose coefficients depend on the PRC 
of the entrained oscillator. A related computation is performed to approximate the 
region in the energy-frequency plane in which entrainment of an oscillator with a given 
PRC by a particular waveform occurs. Such graphs, called Arnold tongues [I9], are 
used to characterize the performance of controls derived using the PRC for entrainment 
of the original oscillator in state space, which is the ultimate purpose of using phase- 
model based control. We conduct this important validation, which is largely lacking in 
the literature, using the Hodgkin-Huxley equations as an example. The results of this 
work can also be viewed as a method for constructing a control to optimally shape the 
Arnold tongue for an entrainment task involving an ensemble of neurons. 

In the following section, we discuss the phase coordinate transformation for a 
nonlinear oscillator and the available numerical methods for computing the PRC. In 
Section [3l we describe how averaging theory is used to study the asymptotic behavior of 
an oscillating system, and use the calculus of variations to derive the minimum energy 
entrainment control for a single oscillator with arbitrary PRC. In Section |U we formulate 
and solve the problem of minimum energy entrainment of oscillator ensembles, for 
which the optimal controls can be synthesized by using an efficient procedure involving 
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Fourier series and Chebyshev polynomials detailed in Appendix A. Throughout the 
paper, important concepts are illustrated graphically by providing examples using the 
Ho dgkin- Huxley model, which is described in Appendix B, and its corresponding PRC. 
We provide computed Arnold tongues in addition to those derived using our theory, 
which verify that the optimal inputs derived here are a significant improvement on 
commonly used entrainment waveforms. In Section [5l we describe several computational 
results that provide further justification for our approach. Finally in Section |6l we 
discuss our conclusions and future extensions of this work. 

2. Phase models 

The phase coordinate transformation is a well-studied model reduction technique that 
is useful for studying oscillating systems characterized by complex nonlinear dynamics, 
and can also be used for system identification when the dynamics are unknown. Consider 
a full state-space model of an oscillating system, described by a smooth ordinary 
differential equation system 

x = f\x,u), x(0) = Xo, (1) 

where x{t) G M" is the state and u{t) G M is a control. Furthermore, we require that 
([I]) has an attractive, non-constant limit cycle 7(t) = 7(t -|- T), satisfying 7 = /(7,0), 
on the periodic orbit F = G M" : y = ^(t) for < t < T} C M". In order to study 
the behavior of this system, we reduce it to a scalar equation 

ij = uj + Z{tlj)u, (2) 

which is called a phase model, where Z is the phase response curve (PRC) and ip{t) 
is the phase associated to the isochron on which x{t) is located. The isochron is the 
manifold in M" on which all points have asymptotic phase ^(t) [39]. The conditions 
for validity and accuracy of this model have been determined |10], and the reduction is 
accomplished through the well-studied process of phase coordinate transformation [H] , 
which is based on Floquet theory 03]. The model is assumed valid for inputs u{t) 
such that the solution x{t,XQ,u) to ([1]) remains within a neighborhood U of F. 

To compute the PRC, the period T = 27r/u; and the limit cycle 7(t) must be 
approximated to a high degree of accuracy. This can be done using a method for 
determining the steady-state response of nonlinear oscillators [B] based on perturbation 
theory [l5] and gradient optimization [16]. The PRC can then be computed by 
integrating the adjoint of the linearization of ([T]) [24], or by using a more efficient and 
numerically stable spectral method developed more recently [1?]. A software package 
called XPPAUT [18] is commonly used by researchers to compute the PRC. The PRC 
of the Ho dgkin- Huxley system with nominal parameters, obtained using a technique 
derived from the method of Malkin [16], is displayed in Figure [TJ 
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Figure 1. Hodgekin- Huxley phase response curve (PRC). The natural period and 
frequency of oscillation arc T « 14.638 and w w 0.429, respectively. 



3. Entrainment of oscillators 



Suppose that one desires to entrain tlie system ([2]) to a new frequency Q using a 
periodic control u{t) = v{Qt) where v is 27r-periodic. We have adopted the weak forcing 
assumption, i.e., v = evi where vi has unit energy and £ << 1, so that given this control 
the state of the original system ([1]) is guaranteed to remain in a neighborhood U of T 
in which the phase model ([2]) remains valid |10]. Now define a slow phase variable by 
(j){t) = ipit) — Qt, and call the difference Au = u — Q between the natural and forcing 
frequencies the frequency detuning. The dynamic equation for the slow phase is 

^ = ip-n = Auj + z{nt + (p)v{nt), (3) 

where is called the phase drift. In order to study the asymptotic behavior of ([3]) it is 
necessary to eliminate the explicit dependence on time on the right hand side, which can 
be accomplished by using formal averaging [18]. Given a periodic forcing with frequency 
fi = 27r/T, we denote the forcing phase 6 = Qt. If V is the set of 27r-periodic functions 
on M, we can define an averaging operator (■) : P — )■ M by 



{x) = - x{e)de. (4) 

Jo 

The weak ergodic theorem for measure-preserving dynamical systems on the torus [17] 
implies that for any periodic function v, the interaction function 



= 1-IJ z{e + <p)v{e)de (5) 
= lim - [ z{nt + ^)v{nt)dt 

T^oo T Jq 

exists as a smooth, 27r-periodic function in V. By the formal averaging theorem [2], the 
system 

^ = Au + A,{ip) + 0{e^) (6) 
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Figure 2. Translated, rescalcd interaction functions of the Hodgkin-Huxley PRC with 
sinusoidal (uo) and square wave {v^) inputs of RMS energy Prms — 0-2, and with 
the optimal waveform u+ for increasing frequency (x) of energy Prms = 0.1301. The 
range of target frequencies 51 to which Vo and w+ can entrain the system are indicated 
by A and B, respectively. Equation ^ has a fixed point for values of up to 2.6% 
greater than ut for each control, but the optimal control requires 35% less RMS energy. 



approximates ([3]) in the sense that there exists a change of variables = + eh{(f, 0) 
that maps solutions of ([3]) to those of ([6]). Therefore the weak forcing assumption v = evi 
with e << 1 allows us to approximate the phase drift equation by 

cp = Auj + A„{ip). (7) 

The averaged equation ([7]) is independent of time, and can be used to study the 
asymptotic behavior of the system ([2]) under periodic forcing. We say that the system 
is entrained by a control u = v{Qt) when the phase drift equation ([7]) satisfies = 0. 
This will eventually occur if there exists a phase ^p^, satisfying Au + A„ = 0. The 
range of frequencies to which the Hodgkin-Huxley system can be entrained by several 
waveforms is illustrated in Figure [2j Because A„((y9) is nontrivial, when the system is 
entrained there exists at least one G [0, 27r) that is an attractive fixed point of ([7]). In 
practical applications, it is desirable to achieve entrainment with a control of minimum 
energy. By defining the phases v?_ = argmin^ At,((/?) and = argmaX;^ A„(93), we can 
formulate minimum energy entrainment of an oscillator as a variational optimization 
problem. The objective function to be minimized is the energy (f^), and entrainment 
can be achieved when uj + Ai,{if^) > Q if Q > u and u + At,(<y5_) < Q if Q < u. This 
inequality is active for the optimal waveform, and hence can be expressed as the equality 
constraint 

Acu + Av{(p+) =0 if n> 

Au; + A^((/7„) =0 if n<u. ^ ^ 

We formulate the problem for Q > u to obtain the minimum energy control w+ using 
the calculus of variations [19]. The derivation of the case where ^2 < a; is similar, and 
results in the symmetric control f„. The constraint ([8]) can be adjoined to the objective 
using a multiplier A, resulting in the cost 
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J[ii\ = {v^)-\{Auj^K{^+)) (9) 



{v^) -x[Auj + — I z{e + ip+)v{e)de 



2^ 



v{9){v{9) - \Z{9 + (^+)) - XAuj]d9. 



Applying the Euler-Lagrange equation, we obtain the necessary condition for an optimal 
solution, which yields a candidate function 

v{9) = ^Z{9 + ^^), (10) 

which we substitute into the constraint (|H]) and solve for the multiplier, A = —2Auj/ (Z^). 
Consequently the minimum energy controls are 

vM = -^,zi9 + cp^) if n>u, 



v.{0) = -^^Z{9 + ^.) if n<oo. 



(11) 



In practice we omit the phase ambiguity (p+ or ip^ in the solution (fTT!) because 
entrainment is asymptotic. 

In addition to deriving the optimal control, we are interested in viewing the Arnold 
tongue, which is a plot of the minimum root mean square (RMS) energy P^(f2) = \/Jv^ 
required for entrainment of the system by a control v to a given target frequency fl. 
This is accomplished by substituting into ([8]) the expression v{9) = Py{Q)vi{9), where 
Vi is a unity energy normalization of v. This results in 

Aid + A,,^{ip+)p^{n) = if n>u, 

Aw + A„,(^_)P^(fi) = if fi<w, ^ ^ 

which in turn yields the Arnold tongue boundary estimate given by 

_Au;/A^,(v9_) if Q<u. ^^^^ 

The boundaries of the theoretical Arnold tongues for the controls f+ and f „, as well 
as computed values of the RMS forcing energy required to entrain the Ho dgkin- Huxley 
system using these controls, are shown in Figure |3l 

In this section, we have shown that the minimum energy periodic control u{t) = v{9) 
that entrains a single oscillator with natural frequency a; to a target frequency f2 is a 
re-scaling of the PRC, where 9 = Qt is the forcing phase. Observe that this control will 
entrain oscillators with natural frequencies between u and Q to the target Q as well. In 
the following section, we derive the minimum energy periodic control that entrains each 
member of a family of oscillators, all of which share the same phase response curve but 
which can have a natural frequency taking any value on a specified interval, to a single 
target frequency Q. 
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Figure 3. Arnold tongues for optimal waveforms for frequency increase (v+) and 
decrease of the Hodgkin-Huxley oscillator. Theoretical boundaries predicted 



by phase reduction theory are shown as lines, and values computed using 3(a) the 



Hodgkin-Huxley PRC and |3(b)| the Hodgkin-Huxley equations are shown as points. 
These values are computed using a line search over the RMS forcing energy, and 3(a) 
closely approximates [3(b)| 



4. Entrainment of an ensemble of oscillators 



We now extend the above approach to derive a single minimum energy periodic control 
signal v{VLt) that guarantees entrainment for each system in the ensemble of oscillators 

= {ifj = u + Z{iIj)u : we(wi,a;2)} (14) 

to a frequency VL. We call the range of frequencies that are entrained by the control v 
the locking range R[v\ = [u-,u+], and when {ui,U2) C -R[f*] we say that the ensemble 
J-" is entrained. This requirement results in the constraints 

Acj„ = Lu^ — n = —Ay{(p^) < wi — = Aui. 

Note that Q G always holds, because an oscillator with natural frequency 

Q is always entrained when forced with the same frequency. The relationship between 
the interaction function the locking range R[v], and the frequency dispersion {ui, U2) 
of the family J-' is illustrated in Figure HI The objective of minimizing control energy 
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Figure 4. The negative of the interaction function (— A„) of the Hodgkin-Huxley 
PRC 4(a) with a square wave (□) input vs of RMS energy P = 0.2 and |4(b)] with the 
optimal waveform for increasing frequency (x) of energy P = 0.1301. The control 
vs achieves a locking range R[vs\ = [w — 0.0112, w -I- 0.0112] that is symmetric about 
the natural frequency w, while the control achieves a non-symmetric locking range 
— [lo — 0.0041, w -f 0.0112]. For target frequencies $7 near a;2, v+ can entrain 
using 35% less energy. We require Aa;2 < Aa;_|- and Awi > Aw- to guarantee 
entrainment of the ensemble. 



{v"^) given the constraints f|T5|l gives rise to the optimization problem 

min J[k] = {v"^), k eV 
s.t. Aw2 + A„(<^_) < 0, 
-Awi - Ay{(^+) < 0. 



(16) 



Uuj2 < ^ (resp. Ui > Q), then problem ( [T6|) is solved by the control f+ where Au = Aui 
(resp. v„ where Au = Aco'2) in (fTT]) . however these are not the only instances where 
that solution is optimal. See, for instance, case (B) in Figure 5(a) Understanding the 
Arnold tongue that characterizes the entrainment of the ensemble J-" by a control v will 
clarify the condition when (ITTj) is optimal. We derive this condition, which depends on 
the ensemble parameters Ui and UJ2 as well as the target frequency Q, and then consider 
the case in which another class of optimal solutions is superior. These two cases are 



illustrated in Figure 5(b) 
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4-0.1. Case I: A re-scaled PRC is optimal for entrainment of the ensemble. To derive 
the conditions when f lTT]) is optimal, we focus on the use of to entrain J-" to a 
frequency Q, G {001^002) when Aa;+ = ^002 > —Aui, noting that the case where 
Au2 < —Aui = — Aci;_ and v+ is used is symmetric. Because U2 is the natural frequency 
in the ensemble farthest from Q, we use Au = Au2- Then the first constraint in (1161) is 
active, yielding —Au2 = Ay{ip_) = Q — u^, so that u+ = 002 is the upper bound on the 
locking range -R[A;], as desired. It remains to determine A„(y9_|_) = Vt — u^, from which 
we obtain the lower bound U- on R[k]. Let us denote A(p = ip^ — and define 

Q{A^) = {Z{e + Aip)Z{e)). (17) 

We can define an inner product (-,■): P x P M by (/, g) = (fg), so that the Cauchy- 
Schwartz inequality yields |Q(Ay9)| < (Z^) = Q{0). Furthermore, the periodicity of Z 
results in Q{Aip) = {Z{e + Aip)Z{e)) = {Z{e)Z{e - Aif)) = Q{-Aip). Combining (P, 
(|TT|) . and ( ITTI) . we can write 

A,(^) = {z{e + v)v4e)) = - (is) 

resulting in 

- Auj+ = A^(v?_) = -Aws, (19) 

as expected. Observe that Ay{ip) is maximized when Q{ip — ip-) is minimized, and hence 
to find A^((y9_|_) it suffices to find the minimum value of Q{A(p). Assume for now that 
< 0, which is true for the Ho dgkin- Huxley PRC, and typical for Type II neurons. 
The practical considerations for finding are discussed in appendix B. It follows that 

A,(^+) = -^Q*, (20) 

and the lower bound of R[k] is = fl — Ay{ip^). If cj_ < ui, then (^1,^2) C R[k], 
hence the control v_ in (II ip . with Au = 002 — ^, is the minimum energy solution to 
problem f|T6|l . and entrains J-" to the frequency Q. 

Now let us define v{d) = Pjj{u)v{6) where v is the unity energy normalization of 
a control v. Substituting this expression into ( !T5|l . we can obtain the minimum RMS 
energy P„_(a;) required to entrain the member of J-" with a natural frequency u to Q 
using the control V-. Because Ay{ip) = Ay{ip)Py{u), this yields 

Au+ + A^_ {<p-)Py_ (u) =0, uj>Q, 

Au^ + A^_{(p+)Py_{uj) =0, uj<n. ^ ' 

Because = Zj-^/TZP), it follows that Aj;_ (y?) = — (p-), and substituting 
this into fl2Tl) and solving for P^, (u) yields 

1 



PvJuj) 



, (cu — il) if f2 < w, 

^^^25(^_(]) if n>uj. 



(22) 
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The boundaries of the Arnold tongues for v„ and f + in fllip are shown in Figure 
5(a) , and examples of oscillator ensembles for which they are optimal are indicated. The 



optimal entrainment problem is fully characterized by the PRC Z and frequency range 
(ci;i,Ci;2) of J-", as well as the target frequency VL. To determine whether the problem is 
optimally solved by ( ITTll . we derive the decision criterion by combining the definition 
A^((^_^) = -Acj. with (EUD and (IIS]) to obtain 



wr 



(23) 



This determines the boundary of the range of VL in relation to (wi, UJ2) when v_ is optimal, 
and the derivation of the optimal range when f+ is optimal is symmetric. Therefore if 



> 



or 



> 



(24) 



then the minimum energy solution to problem ( IT6|) that entrains each member of J-" is 

z{e + ip+) if A002 < -Awi, 

(25) 



AU2 



if 



In practice we omit the phase ambiguity ip^ or (^9+ in solution ( l25l) . The condition 
2M is illustrated as Case I in Figure 5(b) When fl2^ does not hold, then solution fl25l) 



is not optimal for problem ( 1T6|) . as in example (C) in Figure 5(a) which motivates the 
derivation of the optimal solution ( 137|1 below. 



4.0.2. Case II: A difference of shifted PRCs is optimal for entrainment of the ensemble. 



To solve the ( fT6l) when (|25l) is not optimal, as in (D) in Figure 5(a), we adjoin the 



constraints in f lT6|) to the minimum energy objective function using multipliers /i_ and 
giving rise to the cost functional 



J[v] 



(v^) - /i_(Au;2 + - /i+(-Au;i - 

{v')-fi4Auj2 + {z{e + if^)v{e))) 
-fi4-Auji-{z{e + ip+)v{e))) 

^ I \v{eMe) - fi.zie + v^_) + /i+z(e + ^+)] 



yU_Ac<;2 + /i+A(jJi )d6'. 



Solving the Euler-Lagrange equation yields 



(26) 



vie) = --[fi^zie + V.+) - /i„z(^ + (^_)], 



(27) 
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5(a) Arnold tongues of Hodgkin-Huxlcy neuron ensembles for several 

The span of each 



Figure 5 

controls, where the target is the nominal natural frequency. 



horizontal bar represents the range of natural frequencies of a family J- of oscillators 
with the Ho dgkin- Huxley PRC. The vertical location of each bar represents the 
minimum RMS energy required to entrain J- by the indicated waveform. The control 
Vj^ is optimal for ensemble (A), and w_ is optimal for ensemble (B). The waveform u», 
which is optimal for entraining the family (C) to f2 = 5(^1 +^2), achieves entrainment 
with RMS energy 0.26, which is 21% lower than the RMS energy 0.33 required to do 
so with a sine or square wave (D). |5(b)] The appropriate optimal waveform depends on 
the location of O with respect to {uji^uj^). If ^ is in the Case I region then v- (resp. 
v^) is optimal. Otherwise, we use Uh., which depends on wi, W2 ^^nd fi. Recall also the 
assumption < 0. 



which we substitute back into problem f ITB]) to obtain 

{v') = \{{fi^z{e + - f,.z{e + (28) 
= \fil{z') - lfi+f^^{z{9 + ^^)z{e + ^_)) + \^^^z^) 

A,(^+) = {z{e + ^MO)) = -\i^+{z^) + \i^-Q{^v), (29) 
A,(^_) = {z{9 + ^_)v{e)) = \y^-{z^) - ^/i+g(A^). (30) 
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Substituting fl28|) . f l29|) . and f l30|) into problem f|T6|) simplifies the functional optimization 
problem to a nonlinear programming problem in the variables and Q{^ip), 

namely 

min :r[/x_,/i+,Q(A(^)] = + - i;i+/i_g(A^) 

s.t. Aw2 + iyt^-(^'> - |/W+Q(Av9) < 0, (31) 

-Ac^i + \ii+{Z^) - ^fi^QiAy^) < 0. 

When one of the constraints is not active, then either yU,+ = (resp. /i„ = 0), and 
problem ( l26l) is reduced to problem ([9]) with A = /i_ (resp. A = This occurs 

when condition fl2^ holds, and then solution fl25|) is optimal. In the case that condition 
( 124|) does not hold, it follows that both constraints in problem ( 1311) are active. The 
multipliers can be solved for, yielding 



((z2)-g(A^))((z2) + g(Ay,))^ 

2(Aa;ig(A(^) - Aa;2(^')) 

((Z2)-g(Av.))((Z2) + g(A^))' 



(32) 



For these multipliers, the objective in problem fl3T]) is reduced to function of g = Q{A(p) 
given by 

^r^. ^ (Aa;i(Z^) - Aco^Q? + (Ac^^g - Aco^jZ^))' , 

((Z2)-g)2((z2) + g)2 \ / 

2{Aui{Z^) - Auj2Q){AuiQ - Auj2{Z^))Q 



{{Z^) - QY{{Z^) + QY 
Differentiating the cost f l33|) with respect to Q results in 



(33) 



dg " {{z^) - Qmz^) + Qf 



dJ[Q] _ ^ AiOiAio^Q^ - {Aioj + Aiol){Z^)Q + AiOiAio^jZ^y 



Because condition ( l24|) does not hold, it follows that g > Auji{Z^) /Auj2, and hence 
( IMl) is positive for all admissible values of g(A<y9), so that the cost ( 133|) increases when 
g(A<y9) does. Therefore the objective (133|) is minimized when Q{Aip) is, which occurs 
when Q{A(p) = Q* = Q{Alp^). Therefore the problem (l3Ti) is solved when Aip = Aip^ 
and the multipliers are as in fl32|) . Therefore when 

{Z^ Aio^ 

Q. - AC02 - {Z^) ' ^^^^ 

then Case II is in effect, and the minimum energy solution to problem (fT6|) that entrains 
J-" is 

"(''>=((Z^)-e.)((Z^> + <3/<''^''*> 



Optimal Entminment of Neural Oscillator Ensembles 



14 




(a) 




(b) 



Figure 6. 6(a) Arnold tongues of Hodgkin-Huxley phase-model ensembles forced 
by , w+ , a sine wave, and (shown in |6(b)| , where the target $7 is the nominal 
natural frequency. Lines are theoretical values and points are computed using a line 
search over the RMS forcing energy. Note the difference with respect to Figure |31 
Because the Arnold tongue characterizes entrainment of an ensemble of oscillators to 
a fixed target f2 with re-scaled natural frequency uj/^l varying along the a;-axis, its 
contours "drift" right as Prms increasese instead of left, as when the Arnold tongue 
characterizes entrainment of a single oscillator with fixed natural frequency with re- 
scaled target frequency Q,/uj varying along the x-axis. 



The locking range is exactly R[k] = [001,002], satisfying the entrainment constraints ( [T5l) . 
In practice, we omit the phase ambiguity in the solution by using the equivalent control 



({Z''}-Q.){{Z-'}+Q.) 
with energy 

/ 2\ ^ (Acj? + Aool){Z^) - 2A001A002Q, 
^"""^ i{Z^)-Q.)i{Z^) + Q.) ■ ^ ' 

The Arnold tongue for entrainment using v^, can be generated by applying the 
approach used for f_ in ( 12T|) . The advantage of the appropriate optimal control over 
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a sinusoidal or square wave forcing input is illustrated in Figures and 6(a) In the 
special case where VL = ^{uji + ^2), then the solution ( 137|) reduces exactly to the control 
f* that maximizes the locking range R[k] for a fixed control energy, which is achieved 
when Au+ = Au- [STJ |38]. The theoretical contribution presented here provides a 
clarification of the symmetry properties of that particular case. The Arnold tongues 
that characterize entrainment properties of an ensemble J-" are computed and plotted 
for several controls in Figure 6(a), which clearly demonstrates that can entrain an 
ensemble J-" using a lower RMS energy than that required by a sinusoidal waveform. 



5. Robustness of entrainment to parameter uncertainty 

In this section, we provide the results of several numerical simulations that further 
justify our approach to the entrainment of neural ensembles. In particular, we examine 
the effect of parameter variation on the phase response curve and optimal entrainment 
control for an ensemble of neurons. We also provide a visualization of the uncertainty 
in the entrainment properties of a neuron ensemble that arises due to such parameter 
variation. This is done by computing an Arnold tongue distribution, in which the 
minimum RMS energy required for entrainment of the ensemble to a given target 
frequency f2 is a random variable with a probability density on the positive real 
line, instead of a single value, for each u G {ui,U2). We demonstrate that the 
optimal entrainment waveform is minimally sensitive to variation in underlying system 
parameters, that it is always superior to a generic waveform such as a sinusoid or square 
pulse train, and that its amplitude can be appropriately chosen to entrain the neuron 
with the most problematic parameter set in the ensemble. 

The following sensitivity analysis can be performed to examine the effect of 
parameter variation on the PRC of an oscillator. Suppose that pi,...,pd are the 
parameters that characterize the system dynamics, with nominal values ai, . . . , a^. We 
compute the PRC at each corner of a hypercube 

d 

i=l 

where {Pi, 'ji) is a small confidence interval for Oj. The corner points of V are examined 
in particular in order to analyze the aggregate effect of uncertainty in all parameters. 
If the 2'^ curves so obtained are similar to the nominal PRC, then the optimal controls 
derived using the latter will be near optimal for entrainment of an uncertain ensemble. 
Such a robustness property is important in practical neural entrainment applications, 
because biological oscillators exhibit significant variation from any nominal model. Our 
analysis of the sensitivity of the Ho dgkin- Huxley PRC to parameter variation appears 
in Figure [71 in which we have used (/3j,7i) = (0.98ai, 1.02ai) for the d = 7 parameter 
values in the model. For each of the corner points of V we plot the PRC, normalized to 
unity energy, and find that it does not vary significantly from the nominal curve. This 
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(b) 



Figure 7. 7(a) Phase response curves of the Hodgkin- Huxley neuron model and |7(b)] 
the same curves normalized to unity energy. Each of the 2^ grey curves corresponds 
to a PRC obtained when each of the parameters V/va, Vk, Vl, Vng^ Vk^ 5li ^^'^ ^ 
is perturbed by 2% above or below its nominal value. In plot |7(b)| the PRCs of the 
perturbed systems are similar to the black curve, which is the PRC corresponding to 
the nominal parameter set. 



supports our assumption that the optimal entrainment waveform is minimally sensitive 
to variation in underlying system parameters. 

Although minor variations in Hodgkin-Huxley neuron model parameter values 
do not significantly effect the shape of the PRC, they do have a significant effect 
on the entrainment properties of a neuron ensemble by a fixed control waveform v. 
The resulting uncertainty is visualized as an Arnold tongue distribution, which is the 
probability distribution of the minimum RMS control energy required to entrain an 
ensemble of oscillators, with parameter set distributed on a given probability space 
V, as a function of natural frequency u. In practice, we estimate this empirically for 
a hypercube T? with uniform probability measure by uniformly randomly generating 
samples pk G X' of the parameter for k = 1, . . . , iV, for which we compute the natural 
frequency Uk of the perturbed Hodgkin-Huxley model and the minimum RMS control 
energy Pk{v) required to entrain the k^^ model using v. This results in N samples that 
are plotted on the energy-frequency plane, as shown in Figure 8(a), which displays the 
empirical Arnold tongue distribution, using = 500 samples, for a sinusoidal control 
waveform and for the optimal waveform that maximizes the range of entrainment. 
From visual inspection one concludes that the optimal waveform entrains the perturbed 
model using lower power than the sinusoid in the majority of cases. More distinctly. 
Figure 8(b) displays the ratio between the minimum RMS energy levels required to 
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Figure 8. 8(a) Arnold tongue distributions for an ensemble of Hodgkin-Huxley 
neurons with uniform random parameter variation on {j3i,^i) = (.95(3!^, LOSa^), where 
(ai, . . . ^a^) are nominal values for the parameter set Vato, Vk-^ Vl, gNa^ 9k^ 
and c. The target f2 is the nominal natural frequency. There are N — 500 
randomly perturbed parameter sets used to generate empirical distributions (points) 
for entrainment with a sinusoid Vs (o), and with (•), the optimal waveform for 
entrainment to fi = ^(wi +W2)- Solid lines are theoretical Arnold tongue boundaries. 
|8(b)| The ratio Pk{v*) / Pk{vs) plotted as a function of w/O, the natural frequency 
of the neuron with parameter set pk (rescaled by the target frequency) . The optimal 
waveform requires an average of 22% less RMS energy for entrainment than a sinusoid. 



entrain each parameter set pk using the optimal control and a sinusoid, as a function 
of the natural frequency Uk- Not only is the ratio below unity in most cases, with 
an average of 0.78, but there is also a clear trend line at 0.8 in the frequency range 
0.95 < oj < 0.99, with a much lower ratio for natural frequencies near the target 
Vl. This result strongly supports the assertion that the optimal ensemble entrainment 
waveform derived using our method is superior to traditional waveforms such as the 
sinusoid, not only for phase-reduced models, but also for the underlying non-reduced 
dynamical system model. 

In practice, given a confidence region hypercube V for the parameter set P of a 
neuron ensemble, the PRC can be computed at each corner point of V and can be 
used to approximate the corresponding Arnold tongue when the perturbed system is 
entrained by the optimal waveform for the nominal parameter set. In order to assure 
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a/co 

Figure 9. Theoretical Arnold tongues for entrainment of Hodgkin-Huxley neurons 
with parameter sets at the corner points of T) with (l3i,^i) = (0.98ai, 1.02ai), where 
the target is the nominal natural frequency. The Arnold tongues for the nominal 
parameter set (green) and the corner points (grey) are bounded by those of the best- 
and worst- case scenarios. Each is generated using the PRC Zp for the parameter set 
P- 

robust entrainment of the entire ensemble, the RMS energy of v^, should be chosen 
such that it entrains the oscillator with the worst case parameter set, whose theoretical 
Arnold tongue is indicated by the top dashed line in Figure |9l The nominal oscillator is 
entrained with an RMS energy 24% lower than the worst case scenario. This difference 
is very near the average of 22% less RMS energy that the optimal waveform requires 
to entrain an oscillator with parameter uncertainty. It follows that simply by using 
the waveform v^, instead of a square wave or sinusoid one can significantly enhance the 
likelihood that entrainment of an ensemble will be robust to such parameter variation. 

6. Conclusions 

We have developed a method for minimum energy entrainment of oscillator ensembles to 
a desired frequency using weak periodic forcing. Our approach is based on phase model 
reduction and formal averaging theory. We also derive an approximation of the region 
of entrainability in energy-frequency space for an oscillator ensemble. The entrainment 
of phase-reduced Hodgkin-Huxley neurons is considered as an example problem, and 
Arnold tongues are computed to evaluate the effectiveness of our controls. The results 
closely match the theoretical bounds when the weak forcing requirement is fulfilled. The 
optimal waveforms produce a similar result when applied to the original model, which 
suggests that optimal entrainment controls derived using a phase model are optimal 
for the original system, provided the oscillator remains within a neighborhood of its 
limit cycle. We have justified our approach of considering an ensemble with common 
phase response curve and varying frequency by demonstrating robustness of the optimal 
waveforms to variation in model parameters. The results of our simulations suggest 
that stimuli based on inherent dynamical properties of neural oscillators can result in 
significant improvement in energy efficiency and performance over traditional pulses in 
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practical neural engineering applications. This work furthermore provides a basis for 
evaluating the effectiveness of phase reduction techniques for the control of oscillating 
systems with parameter uncertainty. In the future, we will extend the theory of optimal 
entrainment of oscillator ensembles to the case oi n : m entrainment, where m cycles of 
the oscillator occur for every n control cycles. We will also derive controls that most 
rapidly entrain an ensemble of oscillators with uncertain initial state. The approach 
described is of direct interest to researchers in chemistry and neuroscience, and may 
also be applied to vibration control in engineered systems. 

Appendix A. Numerical Issues 

Because the PRC Z{9) and forcing waveform v{9) are both 27r-periodic, we represent 
them using Fourier series, 

^ oo oo 

^(^) = + X] cos(n^) + ^hn sin(n^). (A.l) 

n=l n=l 
oo oo 

v{e) = 2*^0 + ^ Cn cos{n9) + ^ rf„ sin(?T,6'). (A.2) 

n=l n=l 

Using trigonometric angle sum identities and the orthogonality of the Fourier basis, 
we can express the interaction function A„ as 

^ oo 

A^{lp) = {Z{9 + ^)v{9)) = + - ^[a„c„ + bndn] cos{nLp) 

n=l 

^ oo 

+ 9 5Z [^"^" ~ ""^"] sm{nLp) . (A.3) 

^ n=l 

In addition, for any fi,(p2 G [0,27r), periodicity of Z and v ensures that 
{Z{9 + (pi)v{9 + ip2)) = {Z{9 + (fi — ip2)v{9)). Let y = cos{ip), so that we can write 

a c °° 

A^(V2) = Liy) = + - ^KCn + hndn]Tn{y) 

n=l 

^ oo 

+ - ^[hnCn - dndnWl " V^Univ)- (A.4) 
n=l 

where T„ is the n*'^ Chebyshev polynomial of the first kind, and Un is the n^^ Chebyshev 
polynomial of the second kind. If is a minimizer of iv{y) over y G [—1,1], then 
ip^ = arccos(?/*) is a minimizer of Ay{ip) over ip G [— 7r,7r]. Recall that Q{Alp) = 
{Z{9 + Aip)Z{9)) satisfies Q(A^) = Q{-Aip), so that 

2 1 °° 

QiA^) = ^ + 2 + ^n] cos(nA^). (A.5) 

n=l 
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To find and A</?* it is sufficient to minimize a truncation of tfie series in (lA.Sp in tlie 
most numerically expedient way, wfietfier explicitly as Q{^if) = Az{A^p) as in flA.SP 
or by using the change of variables y = cos(A(^) and minimizing iziv) as in ( 1A.4I) . in 
which case the objective function is a polynomial with compact domain [—1,1] C M. 



Appendix B. Hodgkin-Huxley Model 

The Hodgkin-Huxley model describes the propagation of action potentials in neurons, 
specifically the squid giant axon, and is used as a canonical example of neural oscillator 
dynamics. The equations are 

cV = h + I{t)-gj,MV-VNa)m^-9K{V-V,)n^-gL{V-VL) 

m = am{V){l - m) - bm{V)m, 

h = af,{V)il-h)-bhiV)h, 

h = an{V){l - n) - bn(y)n, 



am{V) = 0.1(\/ + 40)/(l-exp(-(y + 40)/10)), (B.l) 

bmiV) = 4exp(-(l^ + 65)/18), 

af,{V) = 0.07exp(-(V + 65)/20), 

bh{V) = l/(l + exp(-(y + 35)/10)), 

aniV) = 0.01(y + 55)/(l -exp(-(F + 55)/10)), 

bn{V) = 0.125 exp(-(V + 65)/80). 

The variable V is the voltage across the axon membrane, m, h, and n are the ion gating 
variables, is a baseline current that induces the oscillation, and I{t) is the control 
input. The units of V are millivolts and the units of time are milliseconds. Nominal 
parameters are V^a = 50 mV, Vk = —77 mV, Vl = —54.4 mV, Qj^^ = 120 mS/cm^, 
= 36 mS/cm^, = 0.3 mS/cm^, lb = 10 fiA/cm^, and c = 1 yuF/cm^, for which the 
period of oscillation is T = 14.63842 ± 10~^ ms. 



References 

[1] S. Strogatz. Nonlinear Dynamics And Chaos: With Applications To Physics, Biology, Chemistry, 
And Engineering. Studies in nonlinearity. Westview Press, 1 edition, 2001. 

[2] F. Hoppensteadt and E. Izhikevich. Weakly connected neural networks. Springer- Vcrlag, New 
Jersey, 1997. 

[3] F.Hanson. Comparative studies of firefly pacemakers. Federation proceedings, S8{S>):2158-216A, 
1978. 

[4] R. Mirollo and S. Strogatz. Synchronization of pulse-coupled biological oscillators. SI AM Journal 

on Applied Mathematics, 50(6):1645-1662, 1990. 
[5] G. Ermentrout and J. Rinzcl. Beyond a pacemaker's entrainment limit: phase walk-through. 

American Journal of Physiology - Regulatory, Integrative and Comparative Physiology, 246(1), 

1984. 

[6] I. Fischer, Y. Liu, and P. Davis. Synchronization of chaotic semiconductor laser dynamics on 
subnanosecond time scales and its potential for chaos communication. Physical Review A, 62, 
2000. 



Optimal Entrainment of Neural Oscillator Ensembles 



21 



[7] I. Blekhman. Synchronization in science and technology. ASME Press translations, New York, 
1988. 

[8] D. Aronson, R. McGehee, I. Kevrekidis, and R. Aris. Entrainment regions for periodically forced 

oscillators. Physical Review A, 33(3):2190-2192, 1986. 
[9] M. Zalalutdinov, K. Aubin, A. Zehnder, R. Hand, H. Craighead, J. Parpia, and B. Houston. 

Frequency entrainment for micromcchanical oscillator. Applied Physics Letters^ 83(16):3281- 

3283, 2003. 

[10] J. D. Bcrkc, M. Okatan, J. Skurski, and H. B. Eichenbaum. Oscillatory entrainment of striatal 

neurons in freely moving rats. Neuron, 43:883-896, 2004. 
[11] A. Sirota, S. Montgomery, S. Fujisawa, Y. Isomura, M. Zugaro, and G. Buzsaki. Entrainment 

of neocortical neurons and gamma oscillations by the hippocampal theta rhythm. Neuron, 

60:683-697, 2008. 

[12] P. A. Tass. A model of desynchronizing deep brain stimulation with a demand-controlled 

coordinated reset of neural subpopulations. Biol. Cybern, 89:81-88, 2003. 
[13] A. M. Kuncel and W. M. Grill. Selection of stimulus parameters for deep brain stimulation. 

Clinical Neurophysiology, 115:2431-2441, 2004. 
[14] E. Izhikevich and Y. Kuramoto. Weakly coupled oscillators. In Encyclopedia of mathematical 

physics. Elsevier, 2006. 
[15] E. Izhikevich. Dynamical Systems in Neuroscience. Neuroscience. MIT Press, 2007. 
[16] I. Malkin. Methods of Poincare and Liapunov in the theory of nonlinear oscillations. Gostexizdat, 

Moscow, 1949. 

[17] I. Kornfeld, S. Fomin, and Y. Sinai. Ergodic theory: Differentiable Dynamical Systems, volume 

245 of Grund. Math. Wissens. Springer- Verlag, 1982. 
[18] Y. Kuramoto. Chemical Oscillations, Waves, and Turbulence. Springer, New York, 1984. 
[19] A. Pikovsky, M. Rosenblum, and J. Kurths. Synchronization: A Universal Concept in Nonlinear 

Science. Cambridge University Press, 2001. 
[20] S. H. Strogatz. From kuramoto to Crawford: exploring the onset of synchronization in populations 

of coupled oscillators. Physica D, 143(l-4):l-20, 2000. 
[21] M. G. Rosenblum, A. S. Pikovsky, and J. Kurths. Phase synchronization of chaotic oscillators. 

Physical Review Letters, 76(11):1804-1807, 1996. 
[22] H. Hong and M. Y. Choi. Synchronization on small-world networks. Physical Review E, 

65:026139(5), 2002. 

[23] H. M. Pinsker. Aplysia bursting neurons as endogenous oscillators, i. phase-response curves for 

pulsed inhibitory synaptic input. J. Neurophysiology, 40(3). 
[24] B. Ermentrout. Type i membranes, phase resetting curves, and synchrony. Neural Computation, 

8(5):979-1001, 1996. 

[25] J. S. Perlmutter and J. W. Wink. Deep brain stimulation. Annu. Rev. Neurosci., 29:229-257, 
2006. 

[26] L. Good. Control of synchronization of brain dynamics leads to control of epileptic seizures in 

rodents. International Journal of Neural Systems, 19(3):173-196, 2009. 
[27] I. Kiss, I. Zhai, and J. Hudson. Emerging coherence in a population of chemical oscillators. 

Science, 296:1676-1678, 2002. 
[28] S. Nakata, K. Miyazaki, S. Izuhara, H. Yamaoka, and D. Tanaka. Arnold tongue of electrochemical 

nonlinear oscillators. Journal of Physical Chemistry A, 113:6876-6879, 2009. 
[29] F. Hoppensteadt and E. Izhikevich. Oscillatory neurocomputers with dynamic connectivity. 

Physical Review Letters, 82(14), 1999. 
[30] J. Hunter and J. Milton. Amplitude and frequency dependence of spike timing: Implications for 

dynamic regulation. Journal of Neurophysiology, 90:387-394, 2003. 
[31] G. UUah. Tracking and control of neuronal hodgkin-huxley dynamics. Physical Review E, 

79:040901(R), 2009. 

[32] A. Nabi and J. Moehlis. Charge-balanced optimal inputs for phase models of spiking neurons. In 



Optimal Entminment of Neural Oscillator Ensembles 



22 



2009 ASME Dynamic Systems and Control Conference, Hollywood, CA, October 2009. 
[33] I. Dasanayakc and J.-S. Li. Optimal design of minimum-power stimuli for phase models of neuron 

oscillators. Physical Review E, 83:061916, 2011. 
[34] I. Dasanayake and J.-S. Li. Charge-balanced minimum-power controls for spiking neuron 

oscillators. IEEE Transactions on Automatic Control (under review). 
[35] A. Hodgkin and A. Huxley. A quantitative description of membrane current and its application 

to conduction and excitation in nerve. The Journal of Physiology, 117(4), 1952. 
[36] J.-S. Li. Control of Inhomogeneous Ensembles. PhD thesis. Harvard University, Cambridge, MA, 

2006. 

[37] T. Harada, H. Tanaka, M. Hankins, and I. Kiss. Optimal waveform for the entrainment of a 

weakly forced oscillator. Physical Review Letters, 105(8), 2010. 
[38] A. Zlotnik and J. Li. Optimal asymptotic entrainment of phase-reduced oscillators. In 2011 ASME 

Dynamic Systems and Control Conference, Arlington, VA, October 2011. 
[39] E. Brown, J. Moehlis, and P. Holmes. On the phase reduction and response dynamics of neural 

oscillator populations. Neural Computation, 16(4):673-715, 2004. 
[40] D. Efimov and T. Raissi. Phase resetting control based on direct phase response curve. In Preprints 

of the 8th IFAC Symposium on Nonlinear Control Systems, pages 332-337, Bologna, September 

2010. 

[41] D. Efimov, P. Sacre, and R. Sepulchre. Controlling the phase of an oscillator: A phase response 
curve approach. In Joint 48th Conference on Decision and Control, pages 7692-7697, December 
2009. 

[42] L. Perko. Differential equations and dynamical systems. Texts in applied mathematics. Springer, 
2 edition, 1990. 

[43] W. Kelley and A. Peterson. The Theory of Differential Equations, Classical and Qualitative. 
Pearson, 2004. 

[44] T. Aprille and T. Trick. A computer algorithm to determine the steady-state response of nonlinear 

oscillators. IEEE Transactions on Circuit Theory, 19(4):354~360, 1972. 
[45] H. Khalil. Nonlinear Systems. Prentice Hall, 3 edition, 2002. 

[46] A. Peressini, F. Sullivan, and J. Uhl. Mathematics of Nonlinear Programming. Springer, 2000. 
[47] W. Govaerts and B. Sautois. Computation of the phase response curve: A direct numerical 

approach. Neural Computation, 18(4):817-847, 2006. 
[48] B. Ermentrout. Simulating, Analyzing, and Animating Dynamical Systems: A Guide to XPPAUT 

for Researchers and Students. SIAM. 
[49] I. Gelfand and S. Fomin. Calculus of Variations. Dover, 2000. 



